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Abstract.  The  complex  response  of  glassy  polymers  to  high  strain-rate  dynamic  loading  necessitates 
accurate  modeling  of  these  events  for  comparison  with  experiments.  The  strain-rate,  temperature  and 
the  strain  softening  behavior  are  significant  and  must  be  considered  for  large  deformations.  Several 
constitutive  relationships  are  discussed  in  terms  of  their  applicability  to  modeling  PMMA  in  gap-tests. 
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INTRODUCTION 

Poly  methyl  methacrylate  (PMMA)  is  often  used 
as  a  window  material  or  as  a  barrier  between  a 
donor  and  acceptor  explosive  in  shock  sensitivity 
experiments  (e.g.  ‘gap-tests’)  [1-6].  However,  the 
response  of  PMMA  (and  glassy  polymers  in 
general)  to  dynamic  loading  is  highly  nonlinear  due 
to  the  competition  of  softening  (temperature  and 
strain)  and  hardening  (strain-rate  and  entropic) 
effects,  which  have  been  investigated  in  [1-10]. 

Strength  models  incorporating  these  effects  evolved 
from  low  strain-rate  processes  (e.g.  extrusion)  in  [7] 
with  strain  rates  in  the  range  of  1 0"4— 1  s'1.  More 
recent  advances  of  this  model  focus  on  high  strain- 
rate  processes  (e.g.  the  Taylor  impact  test)  with  rates 
as  high  as  105  s'1  [8,  9]. 

However,  for  strain  rates  at  or  above  106  s'1  the 
compressibility  of  the  polymer  becomes  significant 
and  the  constitutive  models  must  incorporate  an 
appropriate  equation  of  state.  One  such  constitutive 
model  has  been  developed  [4]  for  modeling  the 
behavior  of  PMMA  in  the  high  pressure  and  strain- 
rate  regime.  It  accounts  for  the  nonlinear  nature  of 
the  shear  strains  and  the  relaxation  present  in  shock 
wave  profiles  but  does  not  incorporate  thermal 


effects.  A  more  recent  constitutive  model  captures 
the  equation  of  state  and  the  prominent  thermal 
effects  derived  from  the  polymer  chain  structure  [8]. 

This  work  discusses  the  first  steps  toward 
modifying  the  Mulliken-Boyce  strength  model  [7] 
for  shock  wave  experiments.  We  also  present  gap 
tests  using  comp-B  as  the  donor  explosive  and 
PMMA  as  the  ‘attenuator’  with  numerical 
calculations  for  comparison. 

The  Mulliken-Boyce  constitutive  model  is 
appropriate  for  glassy  polymers  under  high  strain 
rate  loading  conditions.  It  has  been  implemented  in 
a  Lagrangian  finite  element  code  (EPIC,  2006 
version  [11])  to  capture  the  rate  and  temperature 
effects  in  polymers  under  high  strain  rates  [8,9].  As 
discussed  in  [5,6]  the  strength  of  PMMA  is 
important  to  include  in  gap  test  simulations  in  order 
to  accurately  predict  the  pressure  at  the  interface  of 
the  PMMA  and  receptor  explosive. 

The  implementation  of  the  Mulliken-Boyce 
model  is  depicted  in  Fig.  1,  which  shows  two 
distinct  branches  denoted  A  and  B.  In  branch  A, 
there  are  two  sub-branches  that  correspond  to  the 
glass  (a)  and  viscoelastic  (ft)  transitions  apparent  in 
the  storage  and  loss  moduli  in  dynamic  mechanical 
analysis  (DMA)  [8].  A  modification  to  this  model 
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a 


Figure  1.  Rheological  schematic  of  the  modified 
Mulliken-Boyce  model  [8]. 


is  the  addition  of  a  viscoelastic  dashpot  to  each  sub¬ 
branch  of  ‘A’.  The  model  was  originally 
constructed  to  use  DMA  data  for  the  elastic  moduli, 
but  mesh  dependent  instabilities  occur  in  the 
calculations  due  to  adjacent  elements  having 
different  moduli  and  sound  speeds  [9].  The 
viscoplastic  dashpots  (denoted  ‘v p')  are  somewhat 
of  a  misnomer  since  a  yield  criterion  for  the  material 
is  not  explicitly  defined,  but  the  notation  is 
maintained  for  consistency  with  the  original  model. 

The  kinematics  of  this  model  require  that  the 
total  deformation  in  each  branch  is  equal 

F=Fla  =  F^  =  F£  (1) 


where  F  is  the  deformation  tensor  and  the  subscripts 
denote  the  branch  shown  in  Fig.  1.  The  deformation 
tensor  may  be  multiplicatively  decomposed  into 
elastic  and  thermoplastic  parts,  F  =  Ye¥thp,  where 
¥thp  =  EthW  and  it  is  assumed  that  the  plastic 
deformation  is  an  isochoric  process  (i.e.  detlF^I  =  1 
[4,7-9]).  The  velocity  gradient  in  branch  A  is 

l,=l>f;l'';p(rJ‘  (2} 

where  (*)  denotes  either  a  or  p.  In  Eq.  (2)  the 
thermoplastic  velocity  gradient  is  defined 

I^=D^+W $+0101  (3) 

where  D'/,p  is  the  thermoplastic  rate  of  deformation 
tensor,  a  is  the  coefficient  of  thermal  expansion,  <9  is 
the  temperature  and  I  is  the  identity  tensor.  It  is 
assumed  that  the  plastic  spin  tensor  W///p  =  0.  In  Eq. 
(3)  the  thermoplastic  rate  of  deformation  tensor  is 
coaxial  with  the  branch  deviatoric  stress, 

D%=%WA,+a0A,l  (4) 


where  Np,=<t^,/|cj^J  is  the  normalized  deviatoric 
stress  tensor  and  fpA,  is  the  plastic  shear  strain  rate, 

.  ( AG..  x  \  (5) 


*5.=2  rL.ex] 


sinh  - 

l  kO  s.*+a„ 


i*P) 


In  Eq.  (5),  y*  is  the  initial  plastic  strain  rate 


parameter,  A GA*  is  the  activation  energy  [12],  k  is 


Boltzmann’s  constant,  6  is  the  current  temperature, 
sA*  is  the  temperature  and  rate  dependent  yield 
strength  and  ^  =|(^  \&A  )/2  is  a  measure  of  equivalent 

stress.  Pressure  dependent  yield  is  also  included  in 
Eq.  (5)  with  the  parameter  a  M  >  0  ■ 

The  temperature  and  rate  dependent  strength 
evolves  according  to  the  relation. 


where  h  is  a  strain  softening  parameter  and  sSSiA* 
denotes  the  ‘preferred’  strength  [9].  In  Eq.  (4)  the 
rate  of  temperature  change  is 

(7) 


e=- 


pc, 


^Dj4a  +  Dj4«  +  tr  [**  +  D Afi  J 


where  o  is  the  Cauchy  stress  tensor,  D;  is  the  first 
term  in  Eq.  (4),  p  and  Cp  are  the  density  and  specific 
heat.  The  overbars  in  Eq.  (7)  indicate  that  the  tensor 
resides  in  the  elastically  unloaded  configuration. 
The  viscoelastic  energy  loss  is  accounted  for  in  Eq. 
(7)  with. 


where  rjA*  is  a  viscosity  parameter,  EA*  is  Young’s 
modulus  and  JA  is  the  relaxation  time.  The 
viscoelastic  parameters  for  the  subranches  are 
constants  found  from  DMA  data  by  a  process 
described  in  [9].  The  elastic  constants  are  taken  as 
the  ‘cold’  values  of  the  DMA  data,  however  the 
model  must  incorporate  compressibility  effects 
above  the  Hugoniot  elastic  limit  (0.75  GPa  [2])  for 
comparison  with  experiments. 

In  the  integration  scheme  for  this  constitutive 
model,  the  Jaumann  rate  was  used  to  update  the 
stresses.  This  stress  rate  may  be  expressed  as  a 
function  of  the  rate  of  deformation  tensor, 

(9) 

using  an  isotropic  tangent  stiffness  modulus  CaJM  ■ 


Branch  ‘B’  in  Fig.  1  is  the  polymer  chain 
network  back- stress  describing  entropic  resistance  to 
molecular  alignment.  The  Arruda-Boyce  8-chain 
model  was  used  for  the  back-stress  [13], 


o°) 

B  34*  V4n) 


where  *Jn  is  the  limiting  (or  locking)  value  of 
stretch,  xpcham  =  yjtr(B)/3  *s  a  measure  of  plastic  stretch, 
CR  is  the  rubbery  modulus  £-1(*)  is  the  inverse 
Langevin  function,  and  B'^f\2t2FFT  is  the  deviatoric 


1238 


Downloaded  17  May  2010  to  163.1.15.216.  Redistribution  subject  to  AIP  license  or  copyright;  see  http://proceedings.aip.org/proceedings/cpcr.jsp 


Table  1.  List  of  parameters  for  the  modified  Mulliken- 
Boyce  model  for  PMMA. _ 


Value 

Units 

Value 

Units 

E 

^ Aa 

3386 

MPa 

$ssc/  sQu 

0.73 

- 

&AJ3 

1748 

MPa 

(Xth 

7.4x1  O'5 

1/K 

V 

0.35 

- 

apA- 

0.26 

- 

e0 

298 

K 

fa,  0 

3.2xl013 

1/s 

l Aa 

1979 

s 

f% 0 

3.6xl04 

1/s 

tAp 

1573 

s 

A  Ga 

3.6xl0‘19 

J 

CR 

14 

MPa 

<3=5. 

<1 

6.4x1  O'20 

J 

ha 

300 

MPa 

N 

4.84 

- 

hp 

0 

MPa 

portion  of  the  left-isochoric  Cauchy -Green 
deformation.  Finally,  the  total  stress  is  composed 
from  each  branch  <7=  <TAa  +  +  &b 

The  values  for  the  parameters  were  found  using 
a  genetic  algorithm  in  MatLab®  and  are  listed  in 
Table  1  [14].  The  parameters  are  fit  using 
experimental  data  at  different  temperatures  and 
(constant)  strain  rates.  It  is  assumed  that  with  these 
parameters  the  model  is  sufficient  to  describe  more 
complex  simulations  such  as  the  Taylor  bar  impact 
experiment  shown  in  Fig.  2. 

Gap  tests  have  been  performed  for  model 
validation.  At  this  stage  it  is  useful  to  investigate 
how  well  calculations  compare  to  the  experimental 
data  using  a  cubic  Mie-Gruneisen  equation  of  state 
and  an  elastic-perfectly  plastic  strength  model.  This 
is  an  important  step  for  validation  of  the 
experimental  data  using  parameters  from  well 
documented  sources  [11]  and  serves  as  a  benchmark 
test  for  the  modified  constitutive  model. 

In  each  experiment  a  3”  diameter  PMMA 
cylinder  is  placed  below  a  1”  thick  pad  of  Comp-B 
explosive.  A  streak  camera  was  used  to  record  the 
location  of  the  shock  front  in  time  and  space  by 
rotating  an  image  of  a  flattened  portion,  or  window, 
of  the  PMMA  onto  a  film  strip.  The  setup  was 
backlit  by  a  PBXN-110  pad  (or  Argon  candle).  The 
streak  images  were  processed  using  an  algorithm 
written  in  IGOR  to  extract  the  location  and  time  (or 
x-t  plots)  of  the  shock  front  for  comparison  with 
numerical  calculations. 

Six  PMMA  samples  were  created  for  these 
experiments.  Each  had  a  flat  portion  machined 


Figure  2.  Von  Mises  stress  distribution  in  a  simulated 
Taylor  impact  test  at  250  m/s.  The  image  is  captured  at 
the  point  of  maximum  compression. 

along  the  side  to  allow  light  through  a  slit  into  the 
camera.  One  of  the  samples  (09045)  was  ‘polished’ 
by  flame  while  the  other  five  were  mechanically 
polished.  Also,  a  crack  was  observed  on  the  bottom 
of  one  of  the  samples  (09041).  Cracks  may  act  a 
source  of  enhanced  energy  dissipation  (hot-spots) 
within  the  material  due  to  crack -tip  instabilities  in 
PMMA  [15-17], 

In  Fig.  3  the  streak  records  are  shown.  The 
time  increases  on  the  horizontal  axis  from  left  to 
right  and  the  shock  front,  which  is  the  uppermost 
bright  portion  of  the  curved  section,  passes  through 
the  cylinder  from  the  bottom  to  the  top  of  the  image. 

Simulated  pressure-time  profiles  at  various 
points  in  the  Comp-B  and  PMMA  are  computed  in  a 
simulation  of  the  experimental  setup  using  EPIC06 
[11].  The  cubic  form  of  the  Mie-Gruneisen  equation 
of  state  is  used  for  this  calculation 

P  =  (k,ju  +  k2f/  +  k3ju3 yi-Tju/2)+rEs(l  +  lu)  V1) 
where  Jl-  p/pQ-\ ,  Es  is  the  internal  energy  per  unit 
volume,  kx  =  5.85xl010,  k2  =  3.541xl0n,  k3  = 
1.74xl010,  and  T  =  0.8  is  the  Gruneisen  coefficient. 
The  yield  strength  was  (Jy  =  6.48  MPa. 

Figure  4  shows  the  experimental  results  taken 
from  the  streak  images  compared  to  the  numerical 
calculation.  The  numerical  part  of  the  figure  is  a 
contour  plot  of  the  pressure  with  the  origin  located 
at  the  interface  of  the  comp-B  pad  and  PMMA. 

It  is  clear  that  curves  (3)-(6)  (labeled  from  left 
to  right)  are  tightly  grouped  together  indicating  a 
repeatable  test  with  this  batch  of  PMMA.  Curve  (2) 
corresponds  to  the  ‘flame-polished’  sample,  which 
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Figure  3.  The  streak  image  shows  the  propagation  of  a 
shock  into  PMMA.  Time  increases  from  the  left  to  right 
and  the  propagation  from  bottom  to  top.  The  shock  front 
is  tracked  for  experimental  and  numerical  comparison. 

clearly  exhibits  a  lower  shock  speed.  Curve  (6)  is 
the  sample  with  a  crack  at  the  bottom,  which  did  not 
influence  the  shock  propagation,  which  is  probably 
due  to  its  location.  Curve  (1)  has  the  lowest  speed 
of  the  group.  It  is  unclear  why  this  particular 
sample  exhibited  such  a  low  shock  speed  since  the 
cylinder  was  mechanically  polished  and  free  of 
imperfections.  Overall,  curves  (3)-(6)  match  the 
numerical  calculation  well  until  the  release  wave 
from  the  edge  of  the  cylinder  reach  its  axis. 
Strength  effects  play  an  important  role  in  predicting 
the  time  of  arrival  of  the  release  wave  to  the  center 
axis  of  the  PMMA  and  the  subsequent  pressure 
change  occurring  as  the  shock  front  arrives  at  the 
interface  of  the  acceptor  explosive  in  gap  tests  [6]. 

In  summary,  a  modified  version  of  the  3-D 
Mulliken-Boyce  model  was  implemented  into  EPIC. 
This  model  is  applicable  to  glassy  polymers  in 
general  and  may  be  used  in  the  future  for 
simulations  with  polymers  used  as  binders  in 
explosive  applications.  The  model  lacks 
compressibility  effects,  which  are  essential  in 
modeling  high-pressure  phenomena.  The  cubic 
Mie-Gruneisen  equation  of  state  predicted  the  shock 
speed  well  but  more  sophisticated  strength  models 
are  necessary  to  better  estimate  the  pressure  wave  at 
the  interface  of  the  acceptor  explosive. 
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Figure  4.  Shock  front  profiles  for  six  different  gap  tests. 

The  data  presented  in  curves  (l)-(6)  is  taken  from  an 

algorithm  implemented  in  IGOR. 
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